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In a multi-terminal device the (electronic) heat and charge currents can follow different paths. 

In this paper we introduce and analyse a class of multi-terminal devices where this property is 
pushed to its extreme limits, with charge and heat currents flowing in different reservoirs. After 
introducing the main characteristics of such heat-charge current separation regime we show how 
to realise it in a multi-terminal device with normal and superconducting leads. We demonstrate 
that this regime allows to control independently heat and charge flows and to greatly enhance 
thermoelectric performances at low temperatures. We analyse in details a three-terminal setup 
involving a superconducting lead, a normal lead and a voltage probe. For a generic scattering 
region we show that in the regime of heat-charge current separation both the power factor and the 
figure of merit ZT are highly increased with respect to a standard two-terminal system. These 
results are confirmed for the specific case of a system consisting of three coupled quantum dots. 

PACS numbers: 72.10.-d, 73.23.-b, 74.25.-q, 84.60.-h 


I. INTRODUCTION 

Increasing the efficiency of thermoelectric materials for 
heat-work conversion is one of the main challenges of 
present-days technology 1 In this context the search 
for efficient nanoscale heat engines and refrigerators has 
stimulated a large body of activity, recently reviewed in 
Ref. [71 Progresses in understanding thermoelectricity at 
the nanoscale will also have important applications for 
ultra-sensitive all-electric heat and energy transport de¬ 
tectors, energy transduction, heat rectifiers and refriger¬ 
ators, just to mention a few examples. One of the keys 
to the success in this field is the ability to modulate, con¬ 
trol, and route heat and charge currents, ideally achieving 
their separate controP^ 4 . This is however by no means 
obvious as the charge and (the electronic contribution to) 
the heat are transported by the same carriers. In two- 
terminal systems, for example, within the linear response 
regime, electrical and thermal currents are strictly inter¬ 
related, as manifested by the emergence at low enough 
temperatures of the Wiedemann-Franz law 15 . Indeed, 
when the temperature is the smallest energy scale in 
the system (that is, if the Sommerfeld expansion holds) 
one finds that the ratio A = K/(GT ), involving the 
electrical G and the thermal K conductances at temper¬ 
ature T is universal and it is given by the Lorenz number 
A 0 = tt 2 /3 ( ks/e ) 2 . The fulfillment of the Wiedemann- 
Franz law, one of the triumphs of Sommerfeld’s theory of 
metals, has important consequences in determining the 
efficiency of thermoelectric engines. For a two terminal 
setup the only way to increase the thermoelectric figure 
of merit ZT = (GS 2 / K)T , the dimensionless parameter 
that fully describes the efficiency for thermoelectric con¬ 
version, is by increasing the thermopower ( S ), which is 
however bounded to be small at small X 7 . 



FIG. 1. The heat-charge current separation scheme. A generic 
scattering region is connected to three reservoirs labeled by 
the letters S (superconducting lead), P (voltage probe) and N 
(normal metal lead). In the main text we assume the super¬ 
conducting reservoir (S) to be the reference. In any case, as 
pointed by the arrows, only charge flows inside lead S whereas 
only heat flows inside lead P. 


In order to achieve a separate control of heat and 
charge currents one should therefore consider more com¬ 
plex (multi-terminal) devices. The key issue in this con¬ 
text is to assess to which extent this control can be 
achieved and what are its possible advantages in ther¬ 
moelectric thermal machines. In this paper we push to 
its extreme this type of control and explore a situation 
where heat and charge currents flow in spatially sepa¬ 
rated parts of the system. Namely, we will enforce that 
one of the terminals allows only charge current and an- 
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other one allows only heat current, and name this regime 
as Heat-Charge Current Separation (HCCS). It is worth 
stressing already at this point that the regime of HCCS 
is realised notwithstanding the fact that the same car¬ 
riers are responsible for heat and charge flow. As we 
will discuss in details in the paper, HCCS can be nat¬ 
urally realised by employing superconducting reservoirs. 
Nonetheless, this is not a strict requirement: indeed, in 
principle one could spatially separate heat and charge 
currents in an all-normal multi-terminal device. In this 
case, however, a fine tuning of the parameters character¬ 
ising the thermoelectric transport needs to be performed 
thus making the device not easy to realise experimentally. 

Recent investigations of multi-terminal setups 16 33 
have shown that such devices offer great potentialities 
in terms of efficient thermoelectric conversion. In the 
majority of these works all but two terminals were con¬ 
sidered as probes, i.e. no net flow of energy and charge 
through them was allowed. In other works a purely 
bosonic reservoir was added to the standard two termi¬ 
nals, only exchanging energy with the system. A generic 
three-terminal setup, where all reservoirs are fermionic 
(possibly exchanging both charges and heat with the sys¬ 
tem) , was considered in Ref. [34] where it was shown that 
the third terminal can be used to increase both the ex¬ 
tracted power and the efficiency of a thermal machine. 
So far the impact of superconducting reservoirs on the 
performance of thermoelectric devices has not been con¬ 
sidered^. We start filling this gap by studying HCCS 
and thermoelectric conversion in a three-terminal hybrid 
normal metal-superconducting device. 

The three-terminal device which implements HCCS, 
pictorially shown in Fig. [TJ is composed of a generic con¬ 
ductor connected to a superconducting reservoir (S), a 
normal metal reservoir (N) and a second normal reser¬ 
voir whose chemical potential is set to inhibit the flow 
of electrical current, thus acting as a voltage probe (P). 
This setup, to which we will refer to as SPN, naturally 
realises heat-charge current separation. Indeed, a volt¬ 
age probe exchanges (on average) by definition only heat 
(energy) with the system, whereas the superconductor, 
being a poor heat conductor for temperatures below the 
gap, can exchange only charges. This way, the heat and 
charge currents, flowing together out of the normal metal 
reservoir (N), are split and driven either towards the volt¬ 
age probe (heat), or towards the superconducting reser¬ 
voir (charge). In the linear response regime this setup 
has the advantage of admitting an effective description in 
terms of a 2 x 2 Onsager matrix, a feature which allows 
inter alia a natural way of comparing its performance 
to that of a standard two-terminal configuration. Using 
the scattering approach, we will show on general grounds 
that this separation allows, in the linear response regime 
and at small temperatures, to greatly enhance the per¬ 
formance of a thermal machine, namely increasing both 
the efficiency and the output power by roughly one or¬ 
der of magnitude with respect to a standard two-terminal 
counterpart. The root of this enhanced efficiency can be 


traced back to the possibility to violate in a controlled 
fashion the Wiedemann-Franz law in the heat-charge sep¬ 
aration regime. On more general grounds it is worth to 
stress that the simultaneous presence of superconducting 
and normal terminals, by selectively controlling the heat 
and charge flows through normal and Andreev scattering, 
introduces new degrees of freedom that are worth being 
explored for thermoelectric conversion. 

The paper is organized as follows: in Section [II] we in¬ 
troduce the necessary formalism and define the regime of 
heat-charge current separation. We then show how this 
regime can be attained by having one of the three ter¬ 
minals in the superconducting state. In order to test the 
performance of this thermal engine we perform an exten¬ 
sive analysis in Section m by varying the properties of 
the scattering region connecting the three reservoirs. By 
properly parametrising the scattering matrix we sample 
randomly the scatterer and compare the efficiency of the 
HCCS thermal machine with that of a “conventional” 

and 


two-terminal setup (Sections 
plete our analysis in Section! 


IIIIA 

ncfby discussing the 


of systems consisting of quantum dots (QDs). The reason 
to study these examples in detail is to show that it is pos¬ 
sible to achieve, in experimentally realisable situations, 
those enhanced performances that we found in the first 
part of Section EH Indeed we see that our theoretical 
findings can be tested with current experimental capabil¬ 
ities. Section IV is devoted to the concluding remarks. 
Some technical details related to the scattering formalism 
in the presence of Andreev scattering are summarised in 
the Appendices. 


II. HEAT-CHARGE CURRENT SEPARATION 

Let us consider a system composed of a conductor at¬ 
tached to three leads. Within the linear response regime 
charge and heat currents are governed by the Onsager 
matrix L via the relation 
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where ) represent the charge (heat) current en¬ 

tering the conductor from lead i, with i = (N, P), see 
Fig. [l] We define the biases Xh = A fii/T — (/i — yiq)/T 
and Xf = ATi/T 2 = (T — T^)/T 2 , where yiq and Ti are 
the chemical potential and temperature, respectively, rel¬ 
ative to reservoir i =N,P and having chosen the reservoir 
S as reference with temperature T and chemical poten¬ 
tial fi. Heat and charge currents flowing in lead S can be 
determined from the conservation of particle and energy 
currents. 

As already mentioned, HCCS consists in spatially sep¬ 
arating heat and charge flows. In the example of Fig. [l] 
heat will only flow in lead P while charge will only flow 
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in lead S. In this section we characterise this regime and 
discuss how to implement it. 

On general grounds HCCS can be realised whenever 
two “probe” terminal^ 36 ! are present, one for the voltage 
and one for the temperature. In fact, a voltage probe 
is a terminal whose voltage is adjusted in order for the 
charge current to vanish, while a temperature probe is 
a terminal whose temperature is adjusted in order for 
the heat current to vanish (see App. 0 . Unlike a volt¬ 
age probe, which is implemented simply by opening the 
electric circuit, making a thermal probe would require 
the ability to control and measure heat currents which 
is still very challenging in practice (although important 
advancements in the measurements of heat currents at 
the nanoscale have been recently achieved, see Refs. 37 
and l38|) . A natural way of realising HCCS is to replace 
the thermal probe with a superconducting lead which in¬ 
trinsically suppresses the heat flow for low enough volt¬ 
ages and temperatures. On the contrary a normal metal- 
superconductor interface is an excellent electrical con¬ 
ductor due to the Andreev process that allows to carry 
charge current in the sub-gap regime. In the following we 
will detail the working principles of this implementation. 

Let us consider Eq. 0 and take the superconducting 
reservoir as the reference. Assuming temperatures much 
smaller than the superconducting gap and using the scat¬ 
tering formalism (see App. 0 one can demonstrate that 
the coefficients on the fourth row (column) of the Onsager 
matrix Eq. 0 are the opposite of the corresponding co¬ 
efficients on the second row (column). This implies that 
= —Jp which, at first order in linear response, yields 
J$ = 0 by virtue of the energy conservation. In other 
words, it is an intrinsic property of the hybrid scattering 
matrix to have vanishing heat current in the supercon¬ 
ducting lead. These observations allow us to simplify the 
Onsager system of equations by eliminating the redun¬ 
dant forth row and column, thus reducing it to a 3 by 3 
problem: 


(J c n\ (L n L 12 L 13 \ fX^\ 

\A = \l 21 L 22 L 23 \ \X t , (2) 

\J C P \L 3 1 I/ 32 L 33 \X£ 


where we have introduced X T = Xjj — Xp (or equiva¬ 
lently AT = A Tn — A Tp). Now we impose the voltage 
probe condition Jp = 0 on reservoir P, which yields 


X^ p = 


L 31 X% r + L 32 X T 
£33 


(3) 


By substituting Eq. (|3| into Eq. Q one obtains a two- 
terminal-like Onsager matrix: 


T c 

J N 

J h 

J N 
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For the sake of simplicity, in the following we drop the 
primes for the Onsager coefficients L\ - in Eq. From 
the definitions of the locaP and non-locaP^ transport 


coefficients, one can introduce (local) conductances and 
(non-local) thermopowers described by the following two- 
terminal-like expressions: 
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Importantly, we can express the efficiency for heat to 
work conversion with the standard two-terminal formula 7 

-x*j° n -l^x^-LuX^x* 

n J h N -L 2 - L 22 X t p ' 1 ’ 

One can also define the figure of merit ZT = (GS 2 / K) T 
and the power factor Q = GS 2 . The former gives infor¬ 
mation about the maxim um efficiency and the efficiency 
at maximum powei^® ^(Wmax) = {vc/2) ZT/(ZT + 2), 
r]c = 1 — AT/T being the Carnot efficiency, while the lat¬ 
ter about the maximum power VE ma x = Q{AT) 2 / 4. With 
these formulas the analogy between the SPN system and 
the two-terminal one is complete, allowing us to compare 
their performance. 


III. HSSC IN HYBRID DEVICES 


As we shall show in this section, the heat-charge sep¬ 
aration implemented through the SPN setup allows to 
control G and K separately. This will be at the ori¬ 
gin of the enhancement of both the figure of merit ZT 
and the power factor Q with respect to the two-terminal 
setup. We will use the Landauer-Biittiker scattering for- 
mahsnLLLiyj which is summarised in App. [b] for multi¬ 
terminal hybrid superconducting systems. We begin our 
analysis by considering low temperatures (within the 
Sommerfeld expansion) and studying a well-defined class 
of scattering matrices. Quasiparticle transmission from 
the normal leads into the superconductor is exponentially 
suppressed and thus can b e ign ored. Thus the scattering 
probabilities entering Eq. (B3) just involve reservoirs N 
and P. 

In the following we will express the conductances (elec¬ 
trical and thermal) as well as the thermopower as func¬ 
tions of the parameters characterising the scattering ma¬ 
trix. The aim is to sample this parameter space in order 
to make a statistical analysis of the thermoelectric per¬ 
formance. Assuming a single channel per spin per lead, 
in the Bogoliubov-de Gennes formalism (see App. 0 the 
total scattering matrix Stot is 8 x 8. Supposing that 
there are no spin-mixing terms, it can be written in a 

diagonal block form S to t = ( q <?/ 1 , where the basis is 
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(c t?A r, c t? p, c| ?iV , c| p, c| }N , c| p, c^, c^p), where the op¬ 
erator c a ^ (4-) destroys (creates) an electron with spin 
a in lead i = (TV, S', P). The matrices S and S' are related 
by the particle-hole symmetry relations (see App. [B}, so 
that assigning the elements of S is sufficient to know the 
whole S tot- For sake of simplicity we will consider sym¬ 
metric unitary matrices. A parametrization of such class 
of matrices is given by 


and 


S 2 


f-p 2 (E)e i ^+^+^ y/1 -p 2 (£) 2 e^\ 

V p 2 (E)e~ ie * T’ 


S = 


(9iSi g 2 S 2 \ 

U 5 2 T flsSab 


(9) 


where Si and S 3 are 2 x 2 symmetric unitary matrices, S 2 
is a 2 x 2 unitary matrix, S 3 = SjSJS 2 and < 71 , g 2 , are 

such that the matrix ( J is unitary. For the sake of 

\92 9 s J 

simplicity we assume the latter to be real, i.e. it can be 
g(E) y/l-g(E)*' 
yi - 5(£) 2 - g(E ) 

explicit the dependence on the energy E. Furthermore, 
we parametrize Si and S 2 as 


written as 


where we made 


/- / 9i(£')e i(01+2/3l) sf\ zr fh(E) 2 e^\ 
Pi(E)e~ ie ' ) 


( 10 ) 


assuming the phases in the matrices to be energy in¬ 
dependent. The last simplification that we impose is 
the following relation between the phases: /5 2 + P 2 = 
/3i + 0\ + f. Within this parametrization we assume 
0i, /?i, 0 2 , /3 2 , and y 2 to be real numbers and g(E), 
pi(E) and p 2 (P) to be real functions of energy such that 
0 < g(P), Pi(P), p 2 (P) < 1 for any P. With this nota¬ 
tion and using the expressions for the Onsager coefficients 
in the Landauer-Biittiker approach given in App. [Bj the 
Sommerfeld expansion yields the following transport co¬ 
efficients [see Eqs. 


c — q_ o„('n'i 2 -1_ °v 1 ~f i/w > _ 

9[> -l-p2(0) 2 +fl(0)2(p 1 (0)2 + (3-2p 1 (0)2)p 2 (0)2 + 2(-l + p 1 (0)2)p 2 (0)4)’ 

? g n 2 Tg(0)p2(0) [ - M0)(-1 + MO) 2 ) [(-1 + Pimm + fl(0)M0K(0)] 

-3p 2 (0) 2 + 3g(0) 2 (—1 + pi (0) 2 + (3 - 2 / o 1 (0) 2 )p 2 (0) 2 + 2(-l + pi(0) 2 )p 2 (0) 4 ) 

5 (0)(-l + /9l ( 0 ) 2 )(-l + 2p 2 (0) 2 )^(0)] 

-3/92(0)2 + 3ff(0) 2 (-l + pi (0) 2 + (3 - 2/?! (0) 2 )/? 2 (0) 2 + 2(-l + /9i(0) 2 )/9 2 (0) 4 ) ’ 


K = 


2tt 2 T 
3 


1 + MO) 2 +5(0) 2 (p 2 (0) 2 — 2/9 2 (0) 4 + pi(0) 2 (l — 2p 2 (0) 2 + 2p 2 (0) 4 )) , 
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(13) 

(14) 


where the primed quantities are derivatives with respect 
to energy. After the choices we made, we are left with 
six parameters [namely, pi(0), p[(0), p 2 (0), p 2 (0), g( 0), 
g'(0)] to control G, S and K. We stress that we do not 
impose time reversal symmetry on the scattering matrix 
S, Eq. ([ 9 ]), although our parametrization gives rise to a 
symmetric Onsager matrix.^ 


At this point we would like to draw the attention to 
the fact that the transport coefficients G, S and K are 
independent in a parameter region that is defined by the 
constraints imposed by the unitarity of the scattering 
matrix and from the Sommerfeld expansion. The inde¬ 
pendence of the transport coefficients can be appreciated 
from the way the six parameters, needed to parametrize 
the scattering matrix, enter Eqs. (12)- (14). Indeed, if 
the value of K in Eq. (14) is fixed, the value of G given 
by Eq. (12) is not automatically determined, but instead 
it can be controlled by exploiting the other parameters. 
The same applies to S when G and K are fixed. 


We shall now discuss how the performance of the SPN 
system depends on these parameters. It is indeed im¬ 
portant to verify if: i) HCCS is an advantage for ther¬ 
moelectric conversion, ii) in the regime of HCCS the en¬ 
hancement of the performance is generic or it requires 
additional fine tuning. In order to assess the above is¬ 
sues we will first analyse to which extent heat and charge 
can be controlled independently and then we will study 
the figure of merit as a function of the parameters char¬ 
acterising the scattering matrix. Since we have to deal 
with six free parameters our analysis will be of statistical 
nature. 


A. Control of heat and charge currents 

Our strategy to test our ability to control the currents 
in the three-terminal device is to use, in the same spirit 
as in the Wiedemann-Franz law, the ratio between heat 
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p2 (0) P2 / (0) 


FIG. 2. Plot of the ratio A/Ao and thermopower S for fixed 
K — 10 (SksT)/7r 2 h as a function of P2(0). Here we show that 
using only one parameter [p 2 (0)] we cannot control separately 
the thermopower and the ratio A = K/(GT). The other 
parameters are: pi(0) = 0.8, g'(0) = 0.001(/csT) _1 , p[ (0) = 
Om(kBT)- 1 and p 2 ( 0) = 0.3(fc s T)W 


FIG. 3. Thermopower S from Eq. (13), with K = 
10 (3/csT)/7r 2 /i as a function of p' 2 (0) [in units of ( k B T ) -1 ]. 
Using the additional degrees of freedom provided by the 
derivatives of the parameters [here we use P 2 ( 0 )], we gain 
the control of the thermopower S without affecting the ra¬ 
tio A. The other parameters are: pi(0) — 0.8, g'(0) = 
0.001 (k B T)~\ and pl(0) = 0.1 (ksT)- 1 . 


and electrical conductances. Using Eqs. (12) and (14) 
we can now calculate A = K/(GT ) which can be seen 
by inspection not to be a constant, hence violating the 
Wiedemann-Franz law. Note that both G and K depend 
only on the coefficients pi(0), ^ 2 ( 0 ), and g(0). In Fig. [2] 
we plot, for a fixed value of iF, the dimensionless ratio 
A/Aq and the thermopower S as functions of the pa¬ 
rameter p 2 (0). More precisely, after fixing K we extract 
from Eq. (H4| the parameter g(0) which is a function of 
iF, pi(0), P 21 0) and substitute it into Eq. (12). Moreover, 
we impose the condition that the next order in the Som- 
merfeld expansion is much smaller than the one we take 
into account, restricting the range of admissible values 
of the other parameters (e.g. />2(0) can at most be 0.7). 
For applying this condition we have to specify the val¬ 
ues of the derivatives g'(0), pi(0) and pMO) even though 
they do not appear in Eqs. ( [T2| ) and (Jm|. We assume 
higher order derivatives to be zero for simplicity. The 
plot shows that A is not a constant, but can be controlled 
by properly tuning the parameters of the scattering ma¬ 
trix. Moreover, Fig. [2] shows that S changes by varying 
P 2 (0) for fixed K. The controllability of the transport 
coefficients can be further increased by fixing the values 
of both K and A using the parameters pi(0) and ^ 2 ( 0 ), 
and tuning the derivatives to change the thermopower. 
This is shown in Fig. [3] where A and S are plotted as a 
function of p' 2 ( 0). Notably S spans a quite large interval 
of values, even changing sign. We conclude that the SPN 
system allows independent control of G, K and S. This 
enhanced control is at the basis of the better performance 
that we are going to describe in the next Section. 


B. Thermoelectric performance 


In this section we compare on a statistical ground the 
thermoelectric performance of the SPN system with that 
of a generic two-terminal normal system by randomly 
generating the parameters of the scattering matrices and 
calculating the corresponding power factor Q and figure 
of merit ZT. Within the low temperature limit (Sommer- 




FIG. 4. (left) Probability histogram of the power factor 
Q = GS 2 [in units of k 2 B /h\ for the SPN system (black curve) 
and for the corresponding normal two-terminal system (red 
curve). The maximum value of Q for the SPN setup is about 
0.5 k%/h, while it is about 0.2 k%/h for the two-terminal 
case, (right) Probability histogram of the figure of merit ZT 
for the SPN system (black curve) and for the corresponding 
normal two-terminal case (red curve). The maximum of the 
SPN setup is just above 0.1, while it is about 0.05 for the 
two-terminal system. 

feld expansion), we perform a numerical simulation gen¬ 
erating the parameters of the scattering matrix, of both 
the two-terminal and the SPN systems. Such parameters 
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are picked within a uniform distribution in the allowed 
ranges given by the conditions imposed by the unitarity 
of the scattering matrhP^ and the Sommerfeld expan¬ 
sion. In Fig. [4] we plot the probability of occurrence of a 
certain value of Q (left panel) and ZT (right panel). The 
plot shows that the SPN system (black histograms) has 
better performance than the normal two-terminal system 
(red histograms) for both the power factor Q and the fig¬ 
ure of merit ZT. Indeed, the maximum value of Q for 
the SPN system is about 0.5 k B /h, while it is about 0.2 
k^/h for the two-terminal one. The maximum of ZT for 
the SPN system is just above 0.1, while it is about 0.05 
for the two-terminal system. In Fig. [5] we plot the cor¬ 
relations between Q and ZT for the same random data. 
Each point in the plot corresponds to a particular realiza¬ 
tion of the scattering matrix of the two-terminal or the 
SPN system, for which the power factor and the figure of 
merit are calculated. 



0 - 1 - 1 - 1 - 

0.00 0.02 0.04 0.06 0.08 0.10 0.12 

ZT 


FIG. 5. Correlation between the values of Q and ZT for 
the same data as for Fig.^] each point corresponds to a given 
random realisation. Red (black) points are relative to the two- 
terminal (SPN) setup. The green dashed curve represents the 
bound of Eq. fl5] ), that holds for both the SPN and the two- 
terminal system. The green solid curve, instead, represents 
the bound of Eq. (16), that is the stronger bound given by 
the unitarity on the two-terminal system. 


Fig- 0 shows that the distribution of points presents a 
triangular-like shape for the two-terminal setup. An up¬ 
per bound on the power factor Q is given by the unitarity 
of the scattering matrix. In fact, Q and ZT are related 
by the thermal conductance as Q = (K/T)ZT and, un¬ 
der the Sommerfeld expansion, K is proportional to the 
probability of transmission of an electron from lead N 
to P, which cannot exceed unity. This yields an upper 
bound given by 


with respect to energy imposed by unitarity of the scat¬ 
tering matrhP^. This implies the following expression for 
the maximum of Q 

Qmax = ~y a ° zt ( 1 + ^t\[j 2 ) ’ (16) 

where c is a given energy scale of the order of &gT. Fur¬ 
thermore, for the two-terminal system the power factor 
Q can take all the values between 0 and Q m ax , thus filling 
the red “triangle” of Fig. [5] On the other hand, in the 
case of the SPN system the points are concentrated just 
below the line of the maximum. This is due to the fact 
that the value of K/T, given by Eq. ( p~4| ), cannot take 
all the values between 0 and 2'K 2 k 2 B / (?>h )because of the 
constraints imposed on, and the relations between, the 
parameters appearing in the expression. 

The bound on the maximum value of ZT, instead, is 
given by the conditions on the higher orders of Sommer¬ 
feld expansion, that here we impose to be at least 10 
times smaller that the leading orders for both the two- 
terminal and the SPN system. It is interesting to no¬ 
tice that for the SPN system the points with the highest 
power factor Q are also the points with the highest figure 
of merit ZT: the maximum power automatically gives 
the maximum efficiency! In particular the points with 
the best thermoelectric performance for the SPN system 
roughly correspond to the following values of the scatter¬ 
ing probabilities: normal reflection in lead N, R 0; nor¬ 
mal transmission from lead N to lead P, T — \; Andreev 
reflection in lead N, Ra — A; and Andreev transmission 
from lead N to lead P, Ta — j$- 


C. Coupled QDs in the SPN setup 


We complete our analysis by assessing the thermo¬ 
electric performance of a specific SPN system composed 
of three coupled single-level (non-interacting) QDs (“tri¬ 
dot”) connected to a normal lead (N), a voltage probe (P) 
and a superconducting lead (S), see Fig. [6] The Hamil¬ 
tonian describing the “tridot” reads 


H = 




(17) 


where H e is the Hamiltonian relative to the electrons 
degree of freedom, and is given by 


Q < 
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This is actually true for both the two-terminal and the 
SPN systems. For the two-terminal system a stronger 
bound, which produces the curvature of the upper side 
of the “triangle”, is given by the constraint that must be 
imposed on the derivative of the transmission amplitude 


where q, i = (1,2,3), is the onsite energy of the i —th 
dot, while Tj, {h j} = (1,2,3), is the coupling between 
dot i —th and dot j— th. Notice that Hh = —#* is the 
Hamiltonian relative to the holes, while A = 100 fcgT is 
the superconducting gap. Note that the presence of the S 
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lead is introduced in an effective way, whereby supercon¬ 
ductivity is directly included in the Hamiltonian of the 
“tridot”. The crucial point is that the superconductor 
chemical potential /i is fixed. Furthermore, for simplic¬ 
ity, we have assumed that the superconducting pairing 
for all the QDs is equal as though originating from the 
fact that all QDs are equally coupled to the S lead. In 
Fig. [6j and yp are the coupling energies to the N and 
P lead, respectively. 



FIG. 6. A three coupled QDs (“tridot”) system in the SPN 
setup. 7 at and yp label the coupling to the reservoirs N and P, 
respectively. The three QDs are coupled to a superconducting 
lead (S), with a fixed chemical potential. 


In order to make a statistical analysis we focus on ran¬ 
dom Hamiltonians, for the electron sector i7 e , drawn 
from the Gaussian Orthogonal Ensemble (GOE) and 
from the Gaussian Unitary Ensemble (GUE). The for¬ 
mer describe complex physical systems with time rever¬ 
sal symmetry (TRS), while the latter describe complex 
systems with broken TRfP^HZj The TRS breaking is 
encoded in the complex part of the Hamiltonian (18). 
However this does not imply the spin degeneracy break¬ 
ing, since we can apply a small Aharanov-Bohm flux 
through the plane of the “tridot” (that breaks TRS), 
without a Zeeman component (that would have removed 
the spin degeneracy). In both cases (with and with¬ 
out TRS), the elements of H e are drawn from a Gaus¬ 
sian probability distribution J\f(x,Ax), where x is the 
mean and Ax is the variance. We use distributions 
with different mean for the diagonal elements (QD en¬ 
ergy levels) and for the off diagonal elements (couplings 
between QDs), in order to have a band shift. The for¬ 
mer are drawn from a Gaussian probability distribu¬ 
tion M{k B T , 10 3 k B T) (mean equal to k B T , and variance 
equal to 10 3 £;pT). The latter are drawn from a Gaussian 
probability distribution A7(0, 10 3 k B T). The variance is 
chosen in order to obtain a smooth energy profile for the 
transmission probabilities and to be under the Sommer- 
feld expansion. From the Hamiltonian ( |T7| ) the Green 
function of the system is calculated for fixed values of y^v 
and yp. The 4x4 scattering matrix (there is spin degen¬ 
eracy) is finally calculated from the Green function using 
the Fisher-Lee relation®. Note that the two-terminal 
system, to be compared with the SPN setup, is simply 
described by the Hamiltonian H e . For each random re- 




FIG. 7. Correlation between the power factor Q = GS 2 and 
the figure of merit ZT relative to “tridot” systems for the SPN 
setup (black points) and the two-terminal setup (red points). 
The green curve corresponds to the bound of Eq. ( fl5| , given 
by the unitarity of the scattering matrix and set a maximum 
value for Q as a function of ZT. In panel (a) we show the cor¬ 
relation for the Gaussian Orthogonal Ensemble (GOE), while 
in panel (b) for the Gaussian Unitary Ensemble (GUE). Both 
plots show that for the SPN setup both Q and ZT are one 
order of magnitude larger with respect to the corresponding 
values for the two-terminal system. Moreover it is possible 
to see that the increase of the performance is not due to the 
breaking of the TRS. The plot refers to 10 5 Hamiltonian re¬ 
alizations, taking y/v = yp = y = 10 3 k B T. 


alisation of the Hamiltonian of Eq. © we compute the 
power factor Q and the figure of merit ZT. The results 
are shown in Fig. [7| in both the panels each point in the 
Q — ZT plane represents a single realisation. We notice 
that both Q and ZT are one order of magnitude larger in 
the SPN case with respect to the two-termi nal co nfignra- 
tion for both the GOE Ha milto nian in Fig. 7(a), and the 
GUE Hamiltonian in Fig. |7(b)| This shows a significant 
enhancement, on a statistical ground, in the performance 
of the SPN with respect to the two-terminal system (the 
parameters of the system that realise the maxima of Q 
and ZT for both the two-terminal and the SPN system 




















are given in App. [C|. The distribution of the points is 
similar to that of Fig. [5) Since for the “tridot” model the 
allowed values of ZT (under Sommerfeld expansion) are 
much smaller than the values that we obtained from the 
model in the previous section (see Fig.j5j» , the bound on 
Q of Eq. (16) reduces to that of Eq. (15). Instead, the 
bound on ZT is again set by the Sommerfeld approxi¬ 
mation. Here we notice that for the two-terminal system 
we cannot see the bound given in Eq. (16) because of the 
small values of ZT. As the temperature increases and 
the Sommerfeld expansion loses its validity, we have ob¬ 
served that the performance (for both Q and ZT) of the 
SPN setup gets worse, eventually becoming comparable 
to that of the corresponding two-terminal setup. 


IV. CONCLUSIONS 


In this work we have analyzed the performance of a 
thermal machine which, by involving three reservoirs, al¬ 
lows for the implementation of a spatial separation be¬ 
tween heat and charge currents in linear response. Such 
machine can be naturally realised by connecting a con¬ 
ductor to a superconducting lead, a voltage probe and 
a normal lead (SPN system). Interestingly, the linear- 
response transport equations, written in terms of the On- 
sager matrix, turn out to be formally equal to those of 
a two-terminal conventional system. Using this prop¬ 
erty we have made a comparison between the perfor¬ 
mance of these two thermal machines in terms of the 
power factor Q (that controls the maximum extracted 
power), and the figure of merit ZT (that controls the effi¬ 
ciency at maximum power and the maximum efficiency). 
Within the scattering approach we have described the 
SPN system with a parametrised scattering matrix. We 
have shown that in the low temperature limit (where 
the Sommerfeld expansion holds) the SPN system vio¬ 
lates the Wiedemann-Franz law and allows, to some ex¬ 
tent, an independent control of electrical conductance, 
thermal conductance and thermopower (i.e. of heat and 
charge currents). To assess the consequences of this on 
the thermoelectric performance of the SPN system we 
have made a statistical analysis by taking random values, 
over a uniform distribution, of the parameters contained 
in the scattering matrix. We have thus shown, on statis¬ 
tical grounds, that the SPN system exhibits much larger 
values of Q and ZT with respect to the two-terminal 
counterpart. Further improvements (more than one or¬ 
der of magnitude) of the thermoelectric performance of 
the SPN setup has been confirmed on a specific physi¬ 
cal system composed of three coupled quantum dots. We 
believe that our results can be relevant in the experimen¬ 
tal activity on thermoelectricity of nanoscale structures, 
which are typically conducted at low temperatures. 
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Appendix A: Voltage- and temperature-probe setup 


By imposing the thermal probe condition J = 0 on 
reservoir N and the voltage probe condition Jp = 0 on 
reservoir P, one can solve the second and third rows of 
Eq. Q for the two biases Xj^ and Jp, that will not 
depend directly on the currents J ^ and Jp, and find: 



(l 2 2 L 23 \ 1 (l 21 l 24 \ (x£\ 
\l 32 l 33 ) \l 31 l 34 )\xI)- 


(Al) 


From the first and forth rows of Eq. ([T]) one can define 
an effective two-terminal-like Onsager matrix, with 


Jn 




(A2) 


where the primed Onsager coefficien ts ar e obtained by 
substituting the expressions of Eq. (Al) into Eq. q. 
What has been given above should be considered as a 
definition of HCCS and not a way to implement it. 


Appendix B: Landauer-Buttiker formalism for 
hybrid superconducting systems 

Let us consider a mesoscopic system composed of a 
conductor to which n > 1 leads are attached. Each lead is 
in equilibrium with a fermionic reservoir to which a Fermi 
distribution function is associated, so that a lead is char¬ 
acterized by a temperature and a chemical potential. At 
energy E the i-th lead has Ni(E) open transverse chan¬ 
nels. We allow the possibility to have superconductivity 
in the system and for simplicity we describe it using the 
Bogoliubov-de Gennes (BdG) formalisrrP^, which dou¬ 
bles the degrees of freedom by introducing “hole” states. 
The BdG Hamiltonian is particle-hole symmetric by con¬ 
struction, i.e. it is such that {E[BdG,C} = 0, where C is 
the charge-conjugation operator and the curly parenthe¬ 
ses stand for the anti-commutator. A hole state is the 
charge-conjugate of an electronic state, e.g. if the opera¬ 
tor Ck, a destroys an electron of momentum k and spin cr, 
the operator Cc^C = c\ a destroys a hole of momentum 
k and spin a. For completeness we mention that C is a 
anti-unitary operator hence besides exchanging creation 
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and annihilation operators one must take the complex 
conjugate of the numeric coefficients. 

Particle-hole symmetry implies that the occupation of 
a hole state is the complementary to the occupation of 
an electronic state with opposite energy: 

/-(£) = 1-/+(-£), (Bl) 

where f~ is the distribution function for a hole in lead 
j and is the analogous for electrons. We can then 
write a generalized expression for the Fermi distribution 
function as follows: 


where (ij is the chemical potential of the j —th lead, fi s is 
the chemical potential of the superconductors which we 
take as a reference for the energies, /3j = is the 

inverse temperature of the j -th lead and a is equal to + 
for electrons and — for holes. Assuming coherent trans¬ 
port in the conductor, one can express the charge and en¬ 
ergy currents flowing through the normal leads in terms 
of scattering probabilities using the Landauer-Biittiker 
formalism generalized to include superconductivity: 


ff{E) = 


1 

1 + exp [Pj(E - a(nj - fj, 8 )) 


(B2) 


J? = 


w , p s _ r+oo ' _ r+oo 

= E ( - h) E a / dEP«° & ° (E)/?(E) + - / dENr(E)ft(E), 

n aa0a' ^ ^ 


J“ = 


E(-- h ) E + . 

* aa/3a' JU 


(B3) 


dE(E + a Ms )Nr(-E)ft(E), 


where Jf is the charge current in the i-th lead, J™ is the 
energy current in the i-th lead, e is the electron charge, 
h is the Planck constant and N-* a (E) is the number of 
open channels at energy E for particles of type a and 
spin a. In Eq. (B3) P^ a/3cr (p) is the probability for 
a particle of type /3, spin a' and energy E incoming 
from lead j to be elastically scattered as a particle of 
type a and spin a into the i-th lead. The probabil¬ 
ity of scattering is related to the scattering matrix by 

P^ Pa '{E) = J2 a ,b\S(°aUj,b)(E)\ 2 , where a and b are 
the transverse channels in lead i and j respectively. To 
avoid double-counting that would have been introduced 
by the BdG formalism, the integrals over the energies run 
from 0 to Too instead of starting from — oo. Here zero 


energy corresponds to the Fermi energy of the supercon¬ 
ductors. Due to particle-hole symmetry the probability 
of scattering from lead i to lead j satisfies the relation 

P^’ f}a \E) = Pr j aa ’- f3a \-E). (B4) 

The unitarity of the scattering matrix yields the following 
sum rules: 


E (■ E ) = E p iE°' ( £ ) = n E • 

j,cr',/3 i,cr,a. 

(B5) 

The express ions of Eq . (B3 ) ca n be simplified by substi¬ 
tuting Eqs. (Bl),(B4) and (B5) resulting in 


dE 


(E)Sij Scrcr' 6 0+ - P^\E) 1 ff(E), 


j f=iE r 

E-fe'Eoo 

1 p-\-00 f 

Tt E / dE ( E + N+ a {E)5 ij 5 rja ,5 0+ - P+^ (E) ff(E). 


(B 6 ) 


T U — - 

J% ~ h 


Once the charge and energy currents are determined in 
the normal leads, the sum of the currents in the super¬ 
conducting leads Jsc ~ an d 7 sc = ^7 can 

be calculated exploiting Kirchhoff’s sum rules JE = 0 
and JE = 0 , which are a consequence of charge and 
energy conservation. From the first law of thermodynam¬ 
ics one can also define a heat current j9 in the i-th lead 


as 

= 7 Ji • (B7) 

Let us remark that in a general case there is no sum 
rule for the heat currents. However, one can notice that 
in the superconducting leads the heat current is J- 1 = 
^ Jf, hence the sum of the heat currents over the 
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^sc 


Th _ jh 

~ J SC 


— TU _ Ms TC 

— J SC e J sc 


superconducting leads is 
which can be determined by Kirchhoff’s sum rules on 
charge and energy currents. A hybrid mesoscopic device 
can be thought as a thermal machine, for example for 
heat to work conversion. One can define the efficiency of 
the thermal machine as the ratio between the work W 
extracted from the engine when it absorbs heat Q. By 
convention we assume positive the heat flowing into the 
system, hence in the steady state the definition of the 
efficiency is equivalent to 


rj _W_ 


E jJj _ 

Q E+4* eEtJ? 


(B 8 ) 


where the dot indicates a derivative with respect to time 
and the apex + in the denominator means that the sum 
is restricted to positive heat currents. Time derivative of 
the work W must be positive for the machine to work as 
a heat to work converter, otherwise we are dealing with 
a refrigerator and the definition of 77 is no longer valid. 


J 


Assuming small temperature- and voltage-biases , we 
can expand the Fermi distribution function of Eq. (B2) 
at first order in such quantities: 


f?(E) ~ f{E) + 


dff 


d Xj 


dff 


Y 1 -U 
3 ^ OWM 

(E,T) 


( E,T ) 


v M 

A j > 


(B9) 


dff 


8XJ 


( E,T) 


Edf_ 
T dE' 


dff 




( E,T ) 


,di 

dE' 


where we defined f(E ) = (1 + e fc s T ) _1 . Since the scat¬ 
tering matrix is independent of the biases, we can use 
Eq. ( |B9| ) to linearise the currents of Eqs. (B 6 ) and (BT) 
as follows: 


3 3 


(BIO) 


3 3 

where we defined the quantities 


Gij ~ h 


E f + “ dE [K a (E)S ij 6^ - P+"+"\e) + P^-"'(E)] ( - |1) 

T>H = \ E/ +0 °^ f [N^(E)5 ij 6 aa , - P±° + °\E) - P^-°'(E)\ ( - |£), 

<J<J' °° 

-tZ-.L 

<7 <J 

1 ___ r+o 


M-ij — i 


dEE 


K-a - h 


E 2 r 
dE — 
T 


}( 

df\ 

H 

dE) 

](_ 

df\ 

H 

dEJ 

](- 

df\ 

H 

dE) 


df\ 

H 

dE) 


(Bll) 


Notice that the contributions from energies above the 
gap t o the integrals defining the Onsager coefficients of 
Eqs. (Bll) are exponentially suppressed by the Fermi fac¬ 
tor and will be neglected in the following. Let us notice 
that from the unitarity of the scattering matrix it fol¬ 
lows that the diagonal coefficients are always positive or 
zero. The coefficients ftj, M-ij and /Qj are related 
to the usual Onsager coefficients 7 Lij via the following 
identifications: 


r 


Appendix C: Hamiltonians of the maxima of Q and 
ZT for the “tridot system”. 

In this section we provide the parameters of the elec¬ 
tronic Hamiltonians H e that realise the maximum values 
of Q and ZT for the two-terminal and SPN systems. The 
Hamiltonians that give the maxima for the Gaussian Or¬ 
thogonal Ensemble are 


H: 


2 term, GOE 


1.556 -0.983 —0.046 > 

-0.983 0.109 0.18 

-0.0458 0.18 -0.001, 


7 , (Cl) 


/ An 

A12 

A13 

Al4\ 


( Gn 

TV n 

Q12 

TV i 2 \ 

A21 

A22 

A23 

A24 

— r 

M 11 

TK n 

M .\2 

tk 12 

A31 

A32 

A33 

A34 

— 1 

G21 

TV 21 

G 22 

TV 22 

W 

A42 

A43 

A44/ 


\AI21 

TK 21 

M .22 

TJC 22/ 


(B12) 


'-1.888 -1.074 0.209 
h spn, goe = | 0 74 o 724 -0.834 

0.209 -0.834 0.387 


(C 2 ) 


The Hamiltonians for the maxima for the Gaussian 
Unitary Ensemble are 
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/ 0.729 0.583 + 2.3757 0.273 + 0.115 

H 2 e term - GUE = 0.583 - 2.375 i 0.425 0.062 - 0.226 i 7 , (C3) 

\0.273-0.115i 0.062 + 0.226i 0.039 ) 


^SPN, GUE = 


/ 0.064 -0.923 + 0.484 i -0.16-0.921 A 

-0.923 - 0.484 i 1.213 -0.39 + 0.752 i 7 , 

y—0.16 + 921.13 i —0.39 — 0.752 i 0.481 ) 


(C4) 


where all the energies are expressed in terms of the cou¬ 
pling energy 7 . Note that, although the values of the 
Hamiltonian parameters are much greater than the tem¬ 
perature, the large value of the coupling energy 7 makes 
the linear response coefficients significantly affected by 


such parameters. Indeed, the transmission profiles corre¬ 
sponding to these Hamiltonians turn out to be stretched 
on a (large) scale set by 7 , and hence are almost flat 
within the transport window (that in the linear response 
regime is given by the width of the derivative of the Fermi 
function and is of the order of few fc^T). 
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